Emergent multistability and frustration in phase-repulsive networks of oscillators 
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We study the collective dynamics of oscillator networks with phase-repulsive coupling, considering 
various network sizes and topologies. The notion of link frustration is introduced to characterize 
and quantify the network dynamical states. In opposition to widely studied phase-attractive case, 
the properties of final dynamical states in our model critically depend on the network topology. In 
particular, each network's total frustration value is intimately related to its topology. Moreover, 
phase-repulsive networks in general display multiple final frustration states, whose statistical and 
stability properties are uniquely identifying them. 



I. INTRODUCTION 



Complex systems consist of many individual units 
which generate cooperative functional behavior through 
mutual interactions In the last decade it has been 
realized that complex systems can be elegantly described 
by networks, where nodes represent the functional units, 
and links model their interactions [Ij. The design and 
architectures of networks appearing in nature have been 
extensively studied, revealing a few characteristic classes, 
such as scalefree or modular networks A specific 

emphasis is put on dynamical networks, whose collective 
behavior is a cumulative effect of the individual nodes' 
dynamics and the underlying network topology 0, 0; @1 ■ 
The interplay between network topology and its emer- 
gent dynamics has been widely investigated on various 
examples of empirical and artificially designed networks. 
The coherent dynamics on scalefree networks was found 
to crucially depend on the power-law exponent in the de- 
gree distribution The intra-dependence among the 
dynamical patterns on micro, meso and macro network 
scale was investigated The network constructed by 
coupling its structural evolution to its emergent dynam- 
ics is different from that obtained if these two processes 
are uncoupled Q. Synthetic gene networks are able to 
generate various dynamical regimes in relation to the 
topology of their interactions Q. Collective effects on 
natural networks can be examined by assigning the real 
models as well as formal dynamical systems to the single 
nodes . A recently proposed computational algorithm 
is able to identify the general dynamical patterns induced 
by a given topology independently of the particular dy- 
namical process pTj . 

Models involving repressive or repulsive interactions 
play an important role in the context of dynamical net- 
works. The most popular biological examples are the 
synthetic genetic circuits, in particular toggle switch and 
repressilator [l^] . Genetic circuits consists of a few genes 
that mutually repress each other creating stable oscilla- 
tions of their protein concentrations. Dynamical proper- 
ties of genetic oscillators were extensively studied 0, [ll] . 
Recent works focus on the systems of interacting genetic 
oscillators whose cooperative behavior depends on 

the nature of interactions [l5|, which indicates the ways 
of engineering genetic networks with the desired proper- 



ties. In addition to genes, sparse repulsive cou plin g can 
enhance synchronization in the neural networks jl6j . The 
crucial role of p hase-repulsive interactions was studied 
analytically |17| and confirmed experimentally for oscil- 
lations in neural astrocyte cultures fisj . Repressive in- 
teractions can induce frustration [l^ , which in biological 
systems often generates multistability [l9j . Existence of 
multiple operating regimes is essential for biological sys- 
tems since they provide functional flexibility in respond- 
ing to the stimuli. This has been largely investigated in 
relation to genetic oscillators [TH|, with emphasis on 
biological mechanisms and topological structures lead- 
ing to multistability [20j . The role of multiple dynamical 
regimes was also examined in neuronal interactions, both 
theoretically and experimentally |21| . 

Rhythmic behavior in many natural phenomena can be 
described by a phase variable [22| , which allows the mod- 
eling of complex oscillatory systems and study of the col- 
lective effects such as synchronization (23j . The famous 
Kuramoto model of one-dimensional phase oscillators [l^l 
is widely used not only in theoretical studies [25 Ljbut also 
in modeling specific experimental situations 26|. The 
phase-attractive coupling model was studied in great de- 
tail on a wide range of network sizes and topologies, with 
various distributions of oscillators' frequencies, involv- 
ing different coup ling schemes and time-delayed interac- 
tion [H, 3, H, |23| - |25| . In general, for sufficiently strong 
coupling, the system displays a final synchronized dy- 
namical state, which in the case of identical oscillators 
is always stable full synchronization. Although the time 
scales of the emergence of synchronization may vary [Bj , 
the final network state is in general independent of the 
initial conditions. The time-evolution destroys the infor- 
mation on the network structure, since the most "conve- 
nient" final state always involves synchronization and is 
often completely unrelated to the underlying topology. 

The inherent difference between activatory and re- 
pressory interaction is clearly visible in the phase os- 
cillator models. Phase-repulsive oscillators exhibit alge- 
braic relaxation |27| , in a sharp contrast with the phase- 
attractive case. Despite evolving towards zero mean field, 
arrays of repulsive oscillators display non-trivial dynam- 
ical behaviors such as phase locking and clustering [28j . 
Networks with a given fraction of repulsive links which 
induce dynamical frustration were largely studied [29| . In 
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the context of two-dimensional oscillators, the presence 
of repulsion can improve synchronization [T6| . or even 
generate beam-forming effects that act as a phase array 
antenna [soj . The prescribed synchronization state can 
be achieved through evolutionary network adaptatio n by 
appropriately configuring the repulsive subnetwork |3l|. 
Traveling waves were recently found in a globally coupled 
system with repulsive interactions [s^ . Among many ex- 
perimental scenarios, repulsive oscillators were used to 
model the neuron dynamics with spike ti ming - dependent 
plasticity [s^ and the cultural dynamics 34 1. 

In this paper we consider complex networks of identical 
oscillators with phase-repulsive coupling. Since the oscil- 
lators along each link seek to have the opposite phases, 
wc introduce frustration as a measure of discrepancy with 
this preferred state for each link. Using average link frus- 
tration we characterize the final dynamical states. The 
frustration-topology relationship is systematically ana- 
lyzed employing complex networks of various sizes. In op- 
position to many previous works [l^, H^-fs^ , our network 
model involves only phase-repulsive coupling of uniform 
strength. Considering this simple dynamical model we 
allow for easier study of the interplay between the emer- 
gent dynamics and the underlying topology. As we show, 
contrary to the phase-attractive models, final state of a 
network with phase-repulsive coupling is in general frus- 
trated. The network frustration value is intimately re- 
lated to its topology, and in general increases with its con- 
nectivity. Moreover, repulsive networks in general exhibit 
multiple final dynamical states characterized by different 
values of link frustration. The structure of the frustration 
states directly identifies the network topology, suggesting 
that the repulsive time-evolution preserves much more 
topological information than the attractive one. Repul- 
sive complex networks thus provide a simple model of the 
multistable systems. 

The paper is organized as follows: in the next Section 
we introduce our model and examine its basic properties 
using illustratory examples. In Section IIIII we system- 
atically study all non-directed networks with six nodes, 
analyzing the relationship between topology, connectivity 
and frustration. In Section ITVl we investigate the multi- 
stability of a larger network, considering the transitions 
between frustration states. We discuss our findings and 
conclude in Section IVl 



II. THE MODEL AND BASIC PROPERTIES 

We consider a network consisting of N oscillators 
(nodes) with frequencies LUi. Nodes are connected via 
L non-directed links; — 1 < L < . Dynamical 

state of the oscillator i is described by the phase variable 
(fii € [0, 27r), and its dynamics is given by: 



' J=l,N 



(1) 



where ki is the node's degree {J^i^i = and e is 

the coupling strength. Network's topology is expressed 
through the symmetric adjacency matrix Aij = Aji, with 
value Aij = 1 if nodes i and j are connected, and Aij = 
otherwise. Dynamics starts from a random set of ini- 
tial phases (IP), selected independently for each oscilla- 
tor from (pi{0) G [0, 27r). We consider identical oscillators 
LUi — ij, and take g = sin, thus reducing our system to 
the simple Kuramoto model. Instead of examining the 
model with phase-attractive (positive) coupling e > 0, 
we here focus here on the opposite case involving only 
phase-repulsive (negative) coupling. To this end we fix 
the coupling strength to e = —1. For simplicity we set 
w = 0, i.e. put ourselves in the oscillators' rotating refer- 
ence frame. The equation Eq.(IT]) for our model becomes: 
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A,j sm{ipj - ipi) 



(2) 



The interacting pairs of oscillators are seeking to maxi- 
mize the phase difference between them, i.e. to stretch 
TT apart from each other j27H29|. In the final dynamical 
state, each link will therefore carry the maximal possible 
phase difference, which is preferably tt. However, as we 
show in what follows, due to the complex network topol- 
ogy the phase difference along various links if often less 
than TT, or even zero. 

The global dynamical state of an oscillator ensem- 
ble is usually quantified via the order parameter R = 
W I Sfe 6"^'° I O'^s of its variations for complex net- 
works [11. However, for the purposes of our study, we 
here resort to a different link-based measure of the col- 
lective dynamics. Borrowing the terminology from disor- 
dered systems, we define the frustration fij for each link 
i-j (A„ = 1) as H: 



(3) 



Frustration is related with the impossibility of many in- 
teracting units to simultaneously attain the state of mini- 
mal energy [35j . In our model, a link that stretches to the 
phase difference tt has zero frustration, while for a link 
forced to synchronize [ipj—ifi = 0) the frustration is max- 
imal 2. Frustration measures how "squeezed" is a link: 
it can be pictured as the elastic potential energy con- 
tained in it. We characterize the final (stationary) states 
of dynamical networks Eq.([5]) by assigning a frustration 
value / to each link. To measure the global frustration 
we introduce F as the network average of /: 



^ ~ E/ ^^-J ' 



(4) 



which quantifies how much does the network topology 
allow for links to stretch. F plays the role of non- 
equilibrium potential, since Eq. ^ can be written as [29j : 
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Frustration can be equivalently defined for the phase- 
attractive coupling, with the preferred hnk state having 
zero phase difference. However, since full synchroniza- 
tion is in this case the only final state, all networks will 
trivially have zero F. In contrast, we show here that the 
topology of a phase-repulsive network is reflected in its 
final dynamical state. 

Illustratory examples. We start with the frustra- 
tion on small networks. In our visualization of networks, 
we picture the links using a (color) scale to indicate their 
frustrations fij. The simplest network of two oscillators 
shown in FigU^ is never frustrated: its two nodes always 
attain ipi = and Lp2 = tt, which yields the phase differ- 
ence TT along the link, and hence the frustration f — 0. 
Consider a chain of three oscillators (3-chain) shown in 
FigHt: if the central node has the phase (p2 — 0, nothing 
prevents the other two nodes from having ipi ~ if:} = tt, 
thus again giving the phase difference tt and the frustra- 
tion / = along each link. The situation is however 
different in the case of three-node ring (3- ring): since all 
nodes are now connected, all links can not simultane- 
ously attain the phase difference tt. The stable solution 
is obtained for Lpi = 0,ip2 — ^^fs ~ T"' pl^^-sc 
differences ^ (frustration / = ^) along each link. Due to 
its topology, 3-chain manages to attain a more stretched 
state than 3-ring. 



F=0 

F=0 • • 

> • F=0.5 




F=0.375 F=0.4 F=0.667 

FIG. 1: (color online) Final states of two-node network is 
(a), both three-node networks in (b), and all six four-node 
networks in (c). Total frustration F is reported for each net- 
work. Links are depicted in a (color)scale that indicates their 
frustration fij £ [0,2]. 

In FigjTJ; we show all six four-node networks ordered 
by increasing F. First three of them (top row) achieve 
F ~ due to their specific topologies (we name them 4- 
star, 4-chain and 4-ring, respectively). Note that 4-ring 
in opposition to 3-ring attains F = 0: each diagonal (non 
connected) pair of nodes is synchronized, yielding / = 
along each link. The network with F = 0.375 can be un- 
derstood from the discussion of 3-chain and 3-ring above. 
The fifth network (termed 4-diamond) can be seen as a 4- 



ring with an additional diagonal link. The outside links 
manage to stretch to / = by squeezing the diagonal 
link to / = 2, which gives the total of F = 0.4. Interest- 
ingly, this network achieves the minimal frustration by 
fully squeezing one of its links. The last network is the 
four- node fully connected graph (4-clique) with F ^ ^. 
For different IP this network organizes the values of / 
differently among the links, always achieving the total 
of F = i. The pairs of links that do not share a node 
have the same /-value, and in particular, one such pair 
is always relaxed to / = 0, while the other two divide 
the total of F = ^. The frustration state of 4-clique is 
thus degenerate, with a continuous degeneracy spectrum. 
The total frustration of four-node networks varies with 
both topology and number of links. Each of them has 
a unique way of distributing the frustration among the 
links: while 4-diamond concentrates it into a single link, 
other networks distribute it more uniformly. 

Time-evolution. Contrary to the phase-attractive 
case, time-evolution of the phase-repulsive networks does 
not always exhibit exponential relaxation, and directly 
depends on the network topology [23|. To illustrate this, 
we consider 4-ring and 4-diamond from FigUJ;. For each 
link in each network we examine the behavior of |/(t) — /|, 
where / is the final link frustration value, in addition 
to \F{t) — F\, with F being the final total frustration. 
We show all the curves for a single IP for 4-ring and 4- 
diamond in FigH^ and Figl^Ja, respectively. In the case 
of 4-ring, all four /-values together with the F-value dis- 
play an exponential convergence, similarly to the phase- 
attractive case. In contrast, all five links of 4-diamond ex- 
hibit a power-law convergence with (approximate) slope 
of —1. Interestingly, the convergence of F also shows a 
power-law, but with a steeper slope of —2. 




lO" lO' 10- 10-' 10^ lo' 10- 10-' 

t t 



FIG. 2: (color online) Time-evolution of \f{t)-f\ and \F{t)- 
F\ for 4-ring in (a), and 4-diamond in (b), for a single IP 

These two drastically different convergence regimes re- 
ficct different dynamical processes: while 4-ring quickly 
finds its dynamical equilibrium, the diagonal link of 4- 
diamond resists the phase contraction created by stretch- 
ing of other four links, thus maintaining the system per- 
manently out of equilibrium. These convergence patterns 
are robust to IP, and confirm the earlier findings on the 
relaxation of phase-repulsive oscillators [13]. A similar 
dynamical behavior known as splay states appears in net- 
works of pulse-coupled oscillators [36j . 
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networks 

FIG. 3: (color online) Total frustration values F for all 112 connected six-node networks computed for 10'^ IP. Dashed lines 
divide between groups with different numbers of links L (networks 1-6 have L = 5, networks 7-19 have L = 6 etc.). Within each 
group networks are order by increasing average F (over IP). 11 networks that exhibit muhiple frustration states are marked by 
the vertical (red) lines. 



III. SIX-NODE NETWORKS 

In order to closely examine the relationship between 
the topology and frustration, in this Section we system- 
atically study the phase-repulsive dynamics of all (con- 
nected) six-node networks. There are 112 such networks, 
which wc order by increasing number of links L that 
range from 5 to 15. For each network wc compute the 
total frustration F for lO'^ random IP. Results are re- 
ported in FigI31 where within each group with the same 
L, we order the networks by increasing F (averaged over 
IP), thus constructing a numbering of all six- node net- 
works (numbering serves only to identify the networks). 
Interestingly, there are 11 networks whose final dynami- 
cal state can assume two possible values of F, depending 
on the IP (marked by the vertical lines). The remaining 
101 networks display a unique frustration state, as do 
three-node and four-node networks studied previously. 
The values of F show an overall increase with L, finally 
reaching F = 0.8 for six-node clique. They also exhibit 
large variations within each group with the same L, which 
indicates that the total network frustration strongly de- 
pends on both topology and L. The network group with 
L = 9 links (networks 61-80) exhibits the largest varia- 
tion in frustration values depending on topology; it is also 
the last group where a fully stretched state with = is 
obtainable. Other groups with medium L (from L = 7 to 
L = 10) display the same trend of topology-frustration 
relationship, despite containing different number of net- 
works. Some groups also include many networks with 
similar topologies that all have the same value of F. 

Wc first examine the networks with a unique frustra- 
tion state focusing on the L ~ 9 group. In Fig|4] we 
show the networks 61 and 80 which display minimal and 
maximal _F-values in this group, respectively F — and 
F = 0.611. The network 61 manages to stretch to F = 
by being constructed from symmetrically organized 4- 
stars. Network 80, despite having the same L, has a 
much bigger - it consists of a 3-ring and 4-cliquc, both 
of which have large total frustrations (cf. FiglT|). This 



network exhibits continuously degenerate spectrum of F, 
since it contains the same degeneracy as 4-clique. These 
two opposite examples testify about the flexibility in con- 
taining bigger or smaller frustration within a network, 
realized through variations of its topology. In contrast 



net61,F=0 net 80, F=0.611 




FIG. 4: (color online) Examples of six- node networks dis- 
playing a unique frustration state with values of F indicated. 
Links are marked in (color)scale illustrating their /-values. 

to this, networks 70 and 74 have similar F-values despite 
having rather different topologies. Two networks show 
very different organization of containing the frustration: 
while network 70 distributes it over 8 links, network 74 
confines it into only 2 links, while completely stretching 
the other 7. Network 74 includes two 4-diamonds, which 
in this case display the same frustration pattern as if they 
were isolated. Each topology has its own way of manag- 
ing the frustration, that depends on its particularities 
such as symmetry or modularity. 

Next we study the examples of networks with multiple 
frustration states. In Figl5]we show all 11 of them, visu- 
alized in both states and identified by their numbers as 
described above (cf. Fig|31). For each frustration state we 
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F=0.375 (0.70) net54 F=0.438(0.30) F=0.396 (0.92) net55 F=0.5 (0.08) F=0.432 (0.87) 77 F=0.61 (0.13) 




F=0.545 (0.92) net 102 F=0.727 (0.08) F=0.586 (0.91) net 107 F=0.667 (0.09) 




1 2 



FIG. 5: (color online) All 11 multistable six-node networks visualized in both frustration states and identified by their numbers 
(cf. FiglS}. The F- value is indicated for each state, along with the ratio of IP leading to that state (in parenthesis). (Color)scale 
illustrates links' /-values. 



report the i^- value, along with the fraetion of IP leading 
to it (in parenthesis). There appears to be no specific 
topological property common to all multiple state net- 
works that would distinguish them from the single state 
ones. The simplest multistable network 11 (6- ring) can 
attain the fully stretched state with F = and a squeezed 
state with F = ^. The former is obtained for the phase 
differences tt for all links (equivalently to 4-ring), while 
the latter arises for phase differences ^ along each link 
(equivalently to 3-ring). Both states are stable, but the 
more squeezed one F = ^ occurs for a smaller fraction of 
IP (only 14%). This is to say that both states are stable 
fixed points (sinks) for the dynamical system Eq.([2]) with 
6-ring topology, but the F = state has a bigger basin of 
attraction. Network 35 is a 6-ring with an additional link 
inside: two frustration states differ in distributing the 
frustration between these two subnetworks. Networks 53 
and 54 are topologically similar, and consequently have 
the same F- values occurring for the same fractions of IP. 
Their dynamics is a competition between a 4-ring and 
two 3-rings in escaping the frustration. The F = 0.396 
state of network 55 exhibits a discrete degeneracy: oppo- 
site pairs of links in 4-ring can swap their /-values with- 
out changing F, while the F = 0.5 state shows a con- 
tinuous degeneracy, similar to 4-clique. The remaining 
multistable networks display similar patterns: the differ- 
ence between two states generally lies in the competition 
between two network's structural elements in escaping 



the frustration by attempting to stretch to the maximal 
attainable phase difference. The choice of IP pre-defines 
the final state. The states with lower F- values are usu- 
ally more preferred. The exception to this is network 79: 
its higher F state appears more frequently, despite the 
lower F state involving a uniform distribution of /-values 
over all links. Two F- values are typically close, although 
not always (network 11). Additional stable states are in 
principle possible for some networks, but with extremely 
small basins of attraction, which makes them very diffi- 
cult to observe. Each of the above networks was tested 
for 10^ IP, and no third stable state was found. 

Relaxation of the six-node networks displays the same 
exponential and power-law convergence patterns ob- 
served earlier (Figl2]). Power-law relaxation, testifying 
about the non-equilibrium processes on the network, is 
typically found on the network such as 74 (Fig|3]), which 
concentrate their entire frustration into a few completely 
squeezed links (/ = 2). Interestingly, we revealed various 
power-law slopes for some links in those networks, that 
indicate different squeezing strengths exerted by the rest 
of the network, which is a direct consequence of their 
specific topologies [27| . 

It is instructive to consider the distributions (over IP) 
of initial total frustration F for various networks, since 
they rcfiect their topological symmetries. In FiglS] (top 
panel) we show the distributions of F- values at time t ~ 
for networks 54, 11 and 61 (cf. Figs|3]yS]). The central 



6 



symmetry of distributions for networks 11 and 61 indi- 
cates that each Hnk has the same structural "role" in re- 
lation to the phase-repulsive dynamics. As expected, all 



network 54 network 11 network 61 
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state networks. Phase-repulsive networks almost always 
display many possible final states, whose stability how- 
ever crucially depends on their topology. The appearance 
of multistability can be seen as the persistence of states 
with higher F due to the topological details. 

Smaller networks (of size < 6) do not exhibit multi- 
ple frustration states. On the other hand, multistability 
becomes common as the network size is increased: many 
networks of size N = 7 are multistable, some of them 
possessing three states. As we show in two remaining 
Sections, the number and the organizational complexity 
of frustration states dramatically increases with the net- 
work size and complexity. 



FIG. 6: (color online) Distributions of initial (top panel), and 
final (bottom panel) values of F for many IP, for networks 54, 
11 and 61 (cf. FigsgMSll. 

links of those networks always have the same final values 
of /. On the other hand, the distribution for network 54 
is asymmetric, since not all of its links "see" the network 
in the same way. In the bottom panel of Fig|6]we show 
the distributions of final i^-values, which for networks 54 
and 61 consist of two possible values in a given ratio, 
and for network 61 a unique value F = 0. Note that 
for all networks, the lowest final frustration state is also 
the state with the lowest possible frustration - e.g., for 
network 54, no situation with F smaller than F = 0.375 
is obtainable due to its topology. 

In FigIT] we show the time-evolution of F{t) for 10^ 
IP for network 54. The vertical coordinate indicates the 
fraction of IP having a certain value of F at time t (time 
from t = to t = 20 is considered). We examine the 
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FIG. 7: (color online) Time-evolution of the distribution of F 
for network 54: vertical coordinate indicates the fraction of 
IP (log-scale) with a certain value of F at time t. 10^ IP were 
considered for t G [0, 20]. 

evolution of the initial distribution of F into its two final 
states (cf. Fig l6] left side). Many intermediate unstable 
states with higher F are visited during the evolution be- 
fore settling into one of the final states. For instance, 
a state with F « 0.65 persists for some time, but even- 
tually decays (network 54 was tested for 10^ IP and no 
third stable state was revealed). This cascading dynam- 
ics involving higher frustration states also occurs in single 



IV. FRUSTRATION STATES ON A LARGE 
NETWORK 

In this Section we examine a larger network with more 
frustration states, and study the transitions among them 
occurring by perturbing the dynamics. To this end, we 
construct a network with A^ = 20 nodes as follows: start- 
ing from 3 initially unconnected nodes, we add at each 
step one new node to the existing network. Each new 
node is preferentially attached to two existing nodes, ran- 
domly chosen with probabilities proportional to hi + a, 
where fc^-s are the current node degrees, and a = 1.1. 
The described step is repeated 17 times until the net- 
work size of A^ = 20 is reached, resulting in a 20-node 
network with L = 34 links. Phase-repulsive dynamics 
Eq.(I2]) is implemented as above. 

The dynamics on this preferential attachment grown 
network displays twelve final frustration states. The 
network is visualized in FiglS] in its lowest (left) and 
highest (right) frustration state. We name the states 
Wi,W2, ■ ■ ■Wi2, indexing them by increasing F- value 
termed F{Wi). Each state Wi appears for a certain frac- 
tion of IP which called P{W^). AU 12 values of P{W.,) 
are reported in Figj^^ in relation to the corresponding 
F{Wi). The most preferred state Wi is also the one with 
the lowest F{Wi) = 0.278. The values of P(Wi) overall 
decrease with F{Wi), although the least preferred state is 
Wio with PiWio) < 10-3 and FiWw) = 0.352, while for 
the highest frustration state W12 we find F{Wi2) = 0.364 
and P{Wi2) = 0.006. The states are very unequally 
spaced in F, with F{W4) and ^"(1^5) being nearly the 
same. Each state can be characterized by its specific 
distribution of link frustrations /. As shown in Fig (51 
in state Wi the network stretches most of the periph- 
eral nodes, and confines the entire frustration into the 
links between hubs. In contrast, network in 1^12 stretches 
most of the links around the central hub, while squeezing 
some of the outer links. Similarly, the differences between 
other states typically relate to dividing the frustration be- 
tween the central and peripheral links. In Fig llOl we re- 
port all /-values for all links and all states (numeration of 
links is arbitrary and serves only to discern among them). 
Some links (e.g. 2, 15, 28, 23) exhibit a wide range of 
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F(Wi)=0.278, P(Wj)=0.241 



F(Wj2)=0.364, P(Wj,)=0.006 





FIG. 8: (color online) Lowest (left) and highest (right) frustration state for 20- node network, Wi and W12, visualized with 
(color)scale indicating /. _F-values and P-values are reported. 
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S 0.2 



(a)- 



(b)- 



0.28 0.29 
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FIG. 9: (color online) 12 frustration states Wi,W2, ■ ■ ■W12 
for 20-node network shown through their values of F{Wi). 
(a): fraction of IP P{Wi) leading to each state Wi. (b): the 
fraction of random kicks leading to no change in the state 
PiW, Wi) (cf. Figdlh). 
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FIG. 10: (color online) Link frustration values / for all links 
and all states of 20-node network (arbitrary numeration). Dif- 
ferent symbols are used for different states Wi (legend). 



attainable /-values depending on Wi , that covers the en- 
tire [0,2] interval. Other links (e.g. 7, 11, 13, 24) always 
maintain roughly the same /-value regardless of Wi . The 



former group of links is flexible to different dynamical sit- 
uations, vifhile the latter group is robust to it. Some links 
such as 21 even exhibit two groups of /-values. Some 
pairs of links always have the same /-values for all states 
Wi, which suggests that they have the same dynamical 
role in the network (e.g. and 1, 8 and 20, 12 and 25, 26 
and 30), as also visible in FiglSl The network's response 
to phase-repulsive dynamics involves different dynamical 
roles for different links, realized through a spectrum of 
link frustrations and their flexibility. 

In FigHT] we illustrate the time-evolution of the initial 
distribution of F, as done previously for network 54 (cf. 
FiglZl). The system now visits an even larger number of 
intermediate unstable states with the higher F prior to 
settling in one of the Wi. The speed of this cascading pro- 




FIG. 11: (color online) Time-evolution of initial F distribu- 
tion for 20-node network, as done in Fig[7]for network 54. 10^ 
IP were considered for t G [0, 80]. 

cess seems to decrease with time. The complex topology 
of the underlying network is clearly reflected in the com- 
plexity of general time-evolution. The values of P{Wi) 
are a consequence of starting the dynamics from random 
IP, which yields the initial _F-value much bigger than the 
range of F{Wi). Starting the dynamics from specific IP 
will not influence Wi, but it will change P{Wi). 
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We further compute the Hamming distance H{Wk, Wi) 
between the states defined as: 

HiWk, Wi) = 1 ^ Ay-|/,,(W,.) - f,,{Wi)\ , 

which quantifies the "frustration distance" between any 
two states by averaging the difference in hnk frustrations 
over all links. The symmetric matrix of Hamming dis- 
tances for the 20-node network is shown in FigfT^. The 




FIG. 12: (color online) (a): Hamming distances between frus- 
tration states Wi and Wj for 20-node network, (b): tran- 
sition rates P{Wi — >■ Wj) defined as fractions of random 
kicks yielding a given transition. Note the piecewise constant 
(color)scale. 

states can roughly be divided into four clusters: 

• Wi and W2 have similar distances to all other 
states, and are mutually very close. 

• the same holds for W3, W5 and Wio, which are also 
close to Wi and W2- 

• the cluster of states W7, Wg, Wg, Wn, PF12 are 
higher F states that are mutually relatively close, 
but far from the second, and somewhat close to the 
first cluster. 

• W4 and We are again mutually very close, and also 
close to the first and the third cluster, while far 
from the second. 

The clustering of frustration states is another property of 
phase-repulsive dynamics that reflects the network topo- 
logical details. Note that this classification seems not 
to be directly correlated with the values of F{Wi) and 
P{Wi): for instance, states W4 and W5 are far from each 
other despite having almost the same F-values. 

Below we investigate the transitions between the frus- 
tration states on the 20-node network induced by the 
random perturbations of the network dynamics. To this 
end we modify the Eq.Q by adding the kick term: 

ipi = -— ^ Aij sm{(pj - (fi) + KiSm{ipi + ai)S{t-T) 

' J=l,N 

which acts at time t ^ T by independently perturbing 
the dynamics of each node [25| . For each kick and each 
node, we randomly choose the kicking strength Ki from a 



Gaussian distribution centered at zero with standard de- 
viation 2, and the phase-shifts ai uniformly from [0, 27r). 
The network is prepared at time t = T in state Wi, af- 
ter which the kick is applied. Upon perturbation, the 
network settles into a new state Wj. This procedure is 
repeated 10^ times for each starting state Wi, and the 
transitions Wi — ?> Wj are recorded. We denote with 
P{Wi — ^ Wj) the fraction of random kicks leading from 
the state Wi to the state Wj. The matrix of transitions 
is reported in FiglT2b. where the scale shows the values 
of P{Wi -> Wj). The matrix is of course non-symmetric, 
since it is easier to induce the transitions from a higher 
to a lower F state than vice versa. Similarly, the transi- 
tions from a more preferred into a less preferred state are 
more common than the inverse transitions. In general, 
each state Wi appears to have more and less preferred 
states Wj into which it jumps. The obtained transition 
rates seem to reflect the clustering of the states accord- 
ing to the matrix of Hamming distances shown in Fig |12b : 
the high F states prefer to jump into lower F states that 
belong to their own cluster. For instance, while all states 
jump into Wi and W2 (with various ratios), only some 
states jump into W3, such as members of its cluster W5 
and Wio- Interestingly, among the very few transitions 
occurring from a lower into a higher F state, most start 
from Wi, while all others occur within a given cluster. 

A special role is played by the transitions that do not 
change the frustration state, i.e. Wi — ?> Wi, which are 
the diagonal elements of the matrix P{Wi — >■ Wj) in 
Fig|T2b. They are indicators of the robustness of a given 
state against the perturbations. We show the values of 
P{Wi Wi) in Figinb for comparison with the corre- 
sponding P{W,) in FigHa. The ratios of P{W^ W.,)- 
values only partially reflect the ratios of P(VFi)-values: 
while Wi and W2 arc the most robust states, higher 
F states arc more robust than expected, in particular 
W7,Ws,Wg and Wn. The value of P{Wi) - fraction of 
IP leading to Wi (basin of attraction), can be seen as the 
"width" of the potential hole defining Wi . Similarly, the 
value P{Wi Wi) can be understood as the "depth" of 
the potential hole, as it indicates how strong perturbation 
is needed to jump out of Wi. The comparison of FiglS^ 
and Figinia reveals that depths and widths of the states 
are not completely correlated: low F states are wide and 
relatively deep, while many high F states are only some- 
what shallower despite being much narrower. Note that 
the selection of the kicking strengths is done appropri- 
ately to allow for these properties to be observed. Very 
strong perturbations would erase the memory of start- 
ing state Wi, and all transitions would follow the same 
probabilities as if starting from random IP. 

Finally, we examine the uniqueness of the network frus- 
tration profile (shown in FiglS^) in relation to the net- 
work topology. We implement the link mutation scheme 
as follows: one node of the original network and one of 
its links are chosen at random. The link is then re-wired 
to a different (randomly chosen) node, making sure that 
the network stays connected. The resulting network dif- 
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fers from the original one only in a single link, so it is 
still "topologically close" to it. We compute the statis- 
tics of F(Wi) and P{Wi) for many mutation examples: 
the profile always drastically differs from the original pro- 
file from Figl9^. To illustrate this, we show in FigflBlthc 
original profile (black), together with three examples of 
profiles obtained for networks with a single link mutation. 
The first of them has 20 states, while the second one has 



original network - 
single link mutation 1 - 
single link mutation 2 - 
single link mutation 3 - 



mil. 



0.25 0.26 0.27 0.28 0.29 0.3 0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.4 0.41 0.42 



In the presence of noise our model is expected to ex- 
hibit less frustration states; shallow states such as Wio 
in 20-node network will immediately become unstable. 
With increase of noise strength more states will lose sta- 
bility, finally reaching the point where only a single state 
will remain accessible. This state will thus be the unique 
final dynamical state for a phase-repulsive network. 

For networks with the ring topology the F value is 
directly related to the ring's parity. Recall from FigH] 
that 2-node network and 4-ring have F = 0, while 3- 
ring has F = 0.5. To systematically study this, we show 
in Fig[T3]the F- values of rings as function of number of 
nodes TV. Rings with even TV always have the lowest F 
value F = 0; on the other hand, rings with odd N are 
always frustrated, but their lowest F-value approaches 
zero. With increase of N , rings display a growing num- 



FIG. 13: (color online) 12 frustration states for original 20- 
node network (in black, same as Fig|9^), along with three 
different profiles obtained for three examples of mutated net- 
work (see text for details). 

only 4; the third profile has the most preferred state dif- 
ferent from the lowest F one. It appears that even a 
single link mutation, which only marginally changes the 
topology, yields a dramatic change in the number and 
the properties of the frustration states. This extreme 
sensitivity of multistability to the topology again testi- 
fies about the intricate relationship between them; it ap- 
pears that the frustration profile is in general unique for 
large networks. This can facilitate the reconstruction of 
the phase-repulsive networks from the dynamical data, 
in which context various methods are already in use 



V. DISCUSSION AND CONCLUSIONS 

We studied the collective dynamics of identical Ku- 
ramoto oscillators with phase-repulsive interactions on 
non-directed complex networks. Various network sizes 
and topologies were considered: all 112 connected six- 
node networks were systematically examined, in addition 
to a preferential attachment grown 20-node network. In 
opposition to the phase-attractive case, our model in- 
volves dynamical frustration resulting from the tendency 
of linked oscillator pairs to attain the maximal difference 
of TT between them, which is not always possible due to 
the network's topological complexity. We showed that 
each network has its characteristic total frustration F, 
which largely depends on its size and topology. More- 
over, certain networks display multiple frustration states 
in relation to different initial conditions, which can be 
classified into clusters. Transitions between states also 
refiect topological details and cluster organization. As we 
finally showed, the profile of frustration states appears to 
be a unique "fingerprint" for each network, which is asso- 
ciated with methods of detecting the network structure 
from dynamical data [37| . 



0.5 
Uh 0.4 - 




FIG. 14: (color online) Total frustration F for all rings with 
= 2, . . . 50 nodes computed for 10^ IP. 

ber of multiple frustration states, starting with 6-ring 
(cf. Figl5|). Multiple F- values exhibit a periodic pattern 
with N. The range of attainable F-values shrinks with 
increase of A^, approaching zero at the limit N ^ oo, 
where the ring topology approaches that of a chain. The 
relationship between frustration and parity is associated 
with the methods of detecting network motifs - over- 
represented subnetworks with specific topologies 0, 
Some methods of searching for network rings are already 
in use [11] ■ Since link frustration / contains local net- 
work information, it could be in principle used for motif 
detection. However, this will crucially depend on the way 
motif is embedded in the network - both 4-diamond and 
4-clique (cf. FigH]) contain 3-ring as motif, but its /- 
values arc different from those found on isolated 3-ring. 
One could also seek to generalize the idea of parity in 
the context of networks and find a common topological 
property for all networks with F = 0. 

Another immediate question revolves around the num- 
ber of frustration states in relation to the network con- 
nectivity (number of links L). To examine this we con- 
struct Erdos-Renyi random graphs with A = 40 nodes, 
taking multiples of 10 for L between L = 40 and L = 200 
(L = 40, 50, . . . 200) For each L-value, we construct 

100 different random graph realizations, and record the 
total number of observed states after 200 runs (phase- 
repulsive dynamics is implemented as previously). The 
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results are shown in FiglTS] biggest numbers of states 
(> 50) most often occur on sparse networks around 
L ^ 70. Sparse networks also exhibit the largest range 




40 51) 60 70 KO 'JO 100 110 120 130 140 150 160 170 ISO 190 200 

L 



FIG. 15: (color online) Number of observed F states for each 
realization of Erdos-Renyi random graph with A'' = 40 nodes 
and L links between L = 40 and L = 200. 200 IP were 
considered for each realization. 

of possible number of states depending on the topology, 
and seem always to have no less than 5 states. With in- 
crease of L, networks typically display between 1 and 30 
states, which does not substantially change even for very 
big L. On the other hand, too sparse networks have even 
less states. This result might relate to the sparse con- 
nectivity observed in many biological and technological 
networks 0, 01 j 9'iid emphasize the dynamical proper- 
ties of sparse topologies. This also indicates the optimal 
range of network connectivity for modeling complex mul- 
tistable systems. The 20-node network studied in Section 
HVlis also sparse. 

A further question regards the design of networks with 
minimal or maximal total frustration. We showed in 



Figs [31^4] that a network with a fixed number of links 
may have very different i^-values depending on its topol- 
ogy. It would be interesting to examine the topological 
differences between large networks with fixed L having 
minimal and maximal F. Picturing F as an elastic po- 
tential energy contained in the network, this model may 
indicate the design algorithms for construction of max- 
imally squeezed (or stretched) elastic networks. Similar 
question refers to the networks with minimal or maximal 
number of states, which might be of interest in modeling 
multistable complex systems (cf. Fig fTS]) . 

Future generalizations include networks with non- 
identical oscillators, which are expected to exhibit an 
even wider spectrum of frustration states, including mul- 
tirhythmicity Q. The interaction function g from Eq.([T]) 
was here taken g = sin, although other choices of odd 
g might be interesting. Repulsive dynamics on directed 
and weighted networks is still poorly understood. The 
stability of the fixed points of dynamical system Eq.([2]) 
can also be investigated analytically, using the network 
Laplacian defined as L^j = kiSij — Aij 0,0, HI- Drawing 
conclusions about network multistability by examining 
Lij might allow more detailed and systematic insights. 
In particular, it would be interesting to study the prop- 
erties of the Laplacian eigenvalues in relation to the net- 
work relaxation patterns (cf. Fig 12). 
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